iV-Body Simulations for Coupled Scalar Field Cosmology 
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We describe in detail the general methodology and numerical implementation of consistent N-body 
simulations for coupled scalar field cosmological models, including the background cosmology and 
the generation of initial conditions (with the different couplings to different matter species taken into 
account). We perform fully consistent simulations for a class of coupled scalar field models with an 
inverse power-law potential and negative coupling constant, for which the chameleon mechanism does 
not operate. We find that in such cosmological models the scalar-field potential plays a negligible 
role except in the background expansion, and the fifth force that is produced is proportional to 
gravity in magnitude, justifying the use of a rescaled gravitational constant G in some earlier N- 
body simulations of similar models. We study the effects of the scalar coupling on the nonlinear 
matter power spectra and compare with linear perturbation calculations to investigate where the 
nonlinear model deviates from the linear approximation. For the first time, the algorithm to identify 
gravitationally virialized matter halos is adapted to the scalar field cosmology, and then used to 
measure the mass function and study the properties of virialized halos. We find that the net effect 
of the scalar coupling helps produce more heavy halos in our simulation boxes and suppresses the 
inner (but not the outer) density profile of halos compared with those predicted by lambda-CDM, 
while this suppression weakens as the coupling between the scalar field and dark matter particles 
increases in strength. 

PACS numbers: 04.50.Kd 



I. INTRODUCTION 

The nature of the dark energy [1] driving an apparent 
acceleration of the universe has been a cosmological puz- 
zle for more than a decade. Amongst the models proposed 
to explain it, those incorporating scalar fields are by far 
the most popular, not only because of their mathemat- 
ical simplicity and phenomenological richness, but also 
because the scalar field is a natural ingredient of many 
high-energy physics theories. A scalar field contributes 
a single dynamical degree of freedom which can interact 
indirectly with other matter species through gravity or 
couple directly to matter, producing a fifth force on the 
matter which creates violations of the Weak equivalence 
principle (WEP). This second possibility was introduced 
with the hope that such a coupling could potentially al- 
leviate the coincidence problem of dark energy [2] and 
has since then attracted much attention (see, for exam- 
ple, [3-9] and references therein for some recent work). 

If there is a direct coupling between the scalar field 
and baryons, then the baryonic particles will experience a 
fifth force, which is severely constrained by observations, 
unless there is some special mechanism suppressing the 
fifth-force effects. This is the case in chameleon models, 
where the scalar field (the chameleon) gains mass in high- 
density regions (where observations and experiments are 
performed) and the fifth force effects are confined to un- 
detectably small distances [10-13]. A common approach 
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which avoids such complications is to assume that the 
scalar field couples only to the dark matter, an idea seen 
frequently in models with a coupled dark sector {e.g. [14- 
16]). In this work our scalar field will not be chameleon- 
like as this case has been investigated elsewhere [17-19]. 

The formation of cosmological structure in the pres- 
ence of coupled scalar fields at the linear perturbative 
level has already been studied in great detail. Here, we 
extend these studies into the nonlinear regime. In partic- 
ular, we want to know how structure formation on the 
scales of galaxies and galaxy clusters is modified. As we 
shall see below, there are principally four effects of the 
scalar field coupling, namely through (i) the modification 
of the background expansion rate, (ii) the action of a fifth 
force on dark-matter particles, (iii) the reduced contribu- 
tion of dark matter to the Poisson equation, and (iv) the 
change of the initial condition at early times. Of these 
four, (ii) can actually be subdivided into (a) an essential 
" rescaling" of gravitational constant and (b) a velocity- 
dependent acceleration term, while (iv) is a combined 
consequence of (i, ii, iii) before z ~ 50. We cannot track 
all these effects into the nonlinear regime using linear 
analysis and the tool we will employ to obtain quantita- 
tive predictions is iV-body simulation [20]. 

There has been some earlier work in this area (see 
for example, [21-32] and also [33-38] for related works). 
However, in all those scalar-field simulations the scalar 
field equation is not solved explicitly, but rather it is as- 
sumed either that the fifth force is proportional to gravity 
so that its effect is a simple rescaling of the gravitational 
constant, or that the fifth force takes the Yukawa form 
with exponential cut-off beyond some scale. On the other 
hand, the recent fully consistent simulations performed 



in [17-19] have shown that, at least for the chameleon 
scalar field models, the above simplifying approximations 
arc not good and raises questions about the extent to 
which wc should trust them. It is these concerns that 
motivates this paper, in which we test the accuracy of 
those approximations, and study in more detail the qual- 
itative and quantitative effects of a coupled scalar field 
on the formation and evolution of the nonlinear cosmic 
structure. 

As in [19, 27], we shall consider a universe which at 
late times is dominated by a scalar field and two matter 
species, namely the dark matter, which couples to the 
scalar field in a way prescribed in Sect. II A and the so- 
called "baryons", which are essentially the dark matter 
without scalar-field coupling. 

The paper is organized as follows: in Sect II we list the 
essential equations to be implemented by our TV-body 
simulations and describe briefly their differences from 
standard ACDM. Sect. Ill presents a comprehensive de- 
scription of the methodology used in this paper and its 
implementation in the numerical code. Sect. Ill A intro- 
duces the code, Sect. Ill B lists the physical and simula- 
tion parameters we adopt, Sect. Ill C illustrates how we 
distinguish between baryons and dark matter particles, 
Sect. Ill D summarizes the results for the background ex- 
pansion and linear perturbation evolution, which will be 
referred to subsequently from time to time, and finally 
Sect. HIE and Appendix D explain in detail how to gen- 
erate initial conditions for the A-body simulations which 
take into account the effects of the coupling between dark 
matter and the scalar field. As we aim to set up a general 
framework for A-body simulations of coupled scalar field 
models, we have tried to include all the main ingredi- 
ents in this section. Our numerical results are presented 
in Sect. IV, within which Sect. IV A displays some gen- 
eral results, showing that the approximation of rescaling 
gravitational constant as used in previous literature is a 
very good one for the model studied here (but not nec- 
essarily so for other models!), and Sects. IV B, IV C and 
IV D discuss, respectively, how the coupling modifies the 
nonlinear matter power spectrum, mass function and the 
internal density profiles of halos. We present our conclu- 
sions in V. 

All through the paper a subscript b (d ) denotes a cor- 
responding quantity for baryons (dark matter), unless 
otherwise stated. 



A. The Basic Equations 



The Lagrangian for our coupled scalar field model is 



C 



R 



V>V o¥ > 



■V(<p)-C(<p)CcDu+£a, (1) 



where R is the Ricci scalar, k = 8ttG, with G New- 
ton's constant, ip is the scalar field, V(ip) its potential 
energy, and C(tp) its coupling to the dark matter, which 
is assumed to be cold and described by the Lagrangian 
£cdm; ^s includes all other matter species, including the 
baryons. The contributions from photons and neutrinos 
in the A-body simulations (for late times, z ~ 0(1)) is 
negligible, but should be included when generating the 
matter power-spectrum at redshift z ~ O(50), which de- 
pends on the early evolution of the universe, from which 
the initial conditions for our A-body simulations are ob- 
tained (see Sect. HIE and Appendix D). 

The dark-matter Lagrangian for a point particle with 
(bare) mass mo is 



C 



CDM 



i \ m c-/ / • h 

(y) = — /=<Hy - x )y 9abX%x b , 



(2) 



where y is the general coordinate and xo is the coordinate 
of the centre of the particle. From this equation we derive 
the corresponding energy-momentum tensor: 



rpab 

J CDM — 



m 



9 



5(y-xo)igig. 



(3) 



Also, because gabX^x^ = 9abU a u = 1, where u a is the 
four- velocity of the dark-matter particle centred at xq, 
the Lagrangian can be rewritten as 



n j(, 



£cDM(y) = — y=<5(y - x ), 



(4) 



which will be used below (see Appendix B). 

Eq. (3) is just the energy-momentum tensor for a single 
dark matter particle. For a fluid with many particles, the 
energy-momentum tensor will be 



rpab 
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(5) 



II. THE EQUATIONS 

The equations that will be used in the A-body simula- 
tions have been discussed in detail in Ref. [17, 19]. As we 
are considering a qualitatively different model here, and 
the parameterization and discretization of equations are 
both different, we list these equations and discuss them 
for completeness. 



in which V is a volume that is microscopically large 
but macroscopically small, and we have extended the 3- 
dimensional 5 function to a 4-dimensional one by adding 
a time component. Here, u a is the averaged four-velocity 
of the dark-matter fluid, which is not necessarily the same 
as the four- velocity of the observer, a point which we dis- 
cuss below. 
Using 



rpab 



9 Sg ab 



(G) 



it is straightforward to show that the energy-momentum 
tensor for the scalar field is 



TV ab = ya^yb^ _ g a 



-Vcyvv 



V{<p) 



(7) 



Therefore the total energy-momentum tensor is 



Tab 



\7 a tp\7 b tp - g ab 
+C(^ DM + 






(8) 



where T^ DM = pcDM«a"i), T^ b is the energy-momentum 
tensor for all other matter species including baryons, and 
the Einstein equations are 



^ab — K ^a. 



(9) 



where G a b is the Einstein tensor. Note that because of the 
extra coupling between the scalar field, <p, and the dark 
matter, the energy-momentum tensors for either will not 
be separately conserved, and we have 



V,,T 



CDMafc 






f b C 



CDM 



T 



CDMab 



)Vb<P, 



(10) 

where throughout this paper wc shall use subscript v to 
denote a derivative with respect to tp. However, the total 
energy-momentum tensor is conserved. 

Finally, the scalar field equation of motion (EOM) is 



Dtp- 



dv{<p) dc( v ) 



c 



CDM, 



(11) 



dtp dtp 

where □ = V a V a - Using Eq. (-!-), it can be rewritten as 

dV(tp) dC(tp) 



Dip- 



dtp 



PCDM- 



Op 



= 0. 



(12) 



Eqs. (8, 9, 10, 12) summarize all the physics needed 
for the following analysis. 

We will consider an inverse power-law potential energy 
for the scalar field, 



V(tp) 



A 4 



(y/Ktp) 



(13) 



where a is a dimensionlcss constant and A is a constant 
with dimensions of mass. This potential has also been 
adopted in various background or linear perturbation 
studies of scalar fields (either minimally or non- minimally 
coupled) ; the tracking behaviour its produces makes it a 
good dark energy candidate and for that purpose we shall 
choose a ~ 0(0.1 — 1). Meanwhile, the coupling between 
the scalar field and dark matter particles is chosen as 



C(cp) = exp(7 v / K<^), 



(14) 



where 7 < is another dimensionlcss constant character- 
izing the strength of the coupling. As we shall see below, 



2|7| 2 is roughly the ratio of the magnitudes of the fifth 
force and gravity on the dark matter particles. 

The bare potential Eq. (13) and the coupling function 
Eq. (14) form an effective total potential 



V eff(f) = V(tp) + pcuMC(tp) 



(15) 



for the scalar field tp, just as in [17, 19]. However, both 
V(tp) and C(tp) decrease as tp increases here, there is no fi- 
nite global minimum for V e ff (tp) , and the scalar field will 
always continue rolling down the potential given appro- 
priate initial condition. As a result, although in [17, 19] 
the scalar field almost always resides around the mini- 
mum of V e ff(tp), where it acquires a heavy mass to be- 
come a chameleon [10-13], this does not necessarily hap- 
pen here. Instead, the rolling of the scalar field can be 
quite rapid, introducing interesting new dynamics in both 
background cosmology and perturbation evolution. 



B. The Non-Relativistic Limits 

The A^-body simulation only probes the motion of par- 
ticles at late times, and we are not interested in extreme 
conditions such as black hole formation and evolution, so 
that we can take the non-relativistic limit of the above 
equations as a good approximation. 

As discussed in [19], the existence of the scalar field and 
its couplings to matter particles leads to several changes 
compared with the ACDM paradigm: 

1. The scalar field has its own energy- momentum ten- 
sor, which could change the source term of the Pois- 
son equation because the scalar field, unlike the cos- 
mological constant, can cluster. 

2. The mass of the dark matter particles is effectively 
renormalizcd because of the coupling to the scalar 
field (we say " effectively" because the bare physical 
mass itself is unchanged but it is multiplied by the 
coupling function C(<p) when it appears in gravi- 
tational equations and the scalar field equation of 
motion) . 

3. The scalar field contributes a fifth force on the dark 
matter particles, so they no longer follow geodesies 
determined by gravity only. 

4. As the dark matter and baryonic particles couple 
to the scalar field with different strengths, they will 
influence the scalar field in different ways, feel dif- 
ferent fifth forces, and their subsequent motions will 
differ. 

It therefore becomes clear that the following equations, 
or quantities, in their non-relativistic forms, are needed: 

1. The scalar field equation of motion, to compute the 
value of the scalar field tp at any given time and 
position; 



2. The Poisson equation, to determine the gravita- 
tional potential at any given time and position from 
the local energy density and pressure, which in- 
cludes the contribution from the scalar field (ob- 
tained from ip equation of motion) ; 

3. The total force on the dark matter particles, de- 
termined by the spatial configuration of ip, and like 
gravity it is determined by the spatial configuration 
of the gravitational potential; 

4. The total force on baryonic particles, which receives 
no contribution from the scalar field ip. It is solely 
determined by the spatial configuration of the grav- 
itational potential. 

We shall describe these in turn. For the scalar field 
equation of motion, we denote by ip the background value 
of ip and Sip = <p — (p. Then Eq. (12) can be rewritten as 

Sip + ms'ip + v 2 r p + v v {ip) - v v {ip) 

+Pct>mC }V {p) - PcvmC }V {(p) = 

by removing the background part. Here, V ra is the co- 
variant spatial derivative with respect to the physical 
coordinate r = ax, x is the comoving coordinate, and 
V 2 = V ra Vr- V ra is strictly speaking non-Euclidian as 
the spacetime is not completely flat, but because we are 
working in the weak field limit we approximate it as Eu- 
clidian, that is V 2 = — I d 2 ^ + d 2 + Qj? ) ; the minus sign 

is because our metric convention is (+, — , — , — ). 

In our simulations we also work in the quasi-static 
limit, and assume that the spatial gradients are much 
larger than the time derivatives, |V r </?| 3> |-g^|- There- 
fore, the above equation is simplified to 

c 2 d 2 (aSip) (16) 

= a 3 [V v (ip) - V v ((p) + pcuMC tV {p) - pcbmC, v (ip)} , 

in which d 2 = — V 2 = + (d 2 + d 2 + <9 2 ) is with respect 
to x, with V x = aV r , and we have restored the factor 
c 2 in front of V 2 (the ip here and in the remaining of 
this paper is c~ 2 times the ip in the original Lagrangian 
unless otherwise stated). Note that here V and Pcdm 
both have the dimensions of mass density rather than 
energy density. 

Next consider the Poisson equation, which is obtained 
from the Einstein equation in the weak-field and slow- 
motion limits. Here the metric can be written as 

ds 2 = (1 + 2(f))dt 2 - (1 - 2ip)S ij dr i dr j (17) 

from which we find that the time-time component of the 
Ricci curvature tensor R° = — V 2 cf> . The Einstein equa- 
tion gives 



i? u 



:{ptot + 3ptot) 



(18) 



where ptot and ptot are the total energy density and 
pressure, respectively. The quantity V 2 </> can be ex- 
pressed in terms of the comoving coordinate x as 



1 



$ 1 



-rVi -aax.' 



= ^V 2 $-3^ 
a 3 x a 

where we have defined a new Newtonian potential 

$ = a<j) + -a 2 ax 2 
Thus, 



Vi$ = a J 



V 2 ^- 



(19) 



(20) 



(21) 



-z{Ptot + 3ptot) — -z(Ptot + 3ptot) 



where in the second step we have used Eq. (18) and 
the Raychaudhuri equation, and an overbar labels the 
background value of a quantity. Because the energy- 
momentum tensor for the scalar field is given by Eq. (7), 
it is easy to show that p v + 3p v = 2 [ip 2 — V(ip)] and so 

V 2 $ = -4nGa 3 {pcdmC{p) + p B + 2 [ip 2 - V(ip)] } 



+4irGa 3 {pcdmC(v?) 



Pb 



2[p 2 -V(p)]} 



In this equation ip 2 — ip 2 = 2ipSip + Sip ^C (V r ip) 2 in 
the quasi-static limit, and so could be dropped safely. 
Therefore, we have 



<9 2 $ = AttGo 3 [pcumC{p) - pcumC(p)} 

+AirGa 3 [p B - p B ] - SirGa 3 [V{ip) - V{p)}{22) 

Finally, for the equations of motion of the dark matter 
particles, consider Eq. (10). Using Eqs. (3, 4), this can 
be reduced to 



bc^O^O 



5 



C(p) 



V b ip. (23) 



In this equation, the left-hand side is the conventional 
geodesic equation of general relativity, and the right-hand 
side is the new fifth force contribution from the coupling 
to the scalar field. Note that g ab - u a u b = h ab is the pro- 
jection tensor that projects any 4-tensor into the 3-space 



perpendicular to u a , so (g 



ah 



vru 



') V a = V 6 is the spa- 



tial derivative in the 3-space of the observer and perpen- 
dicular to u a . Therefore, the right-hand side of Eq. 



c(m\ v a <l = V a logC(v?), is perpendicular to u a , which 
shows that the energy density of dark matter is conserved 
and only the particle trajectory is modified [17]. As the 
'observers' are just the individual dark matter particles 
[cf. Eq. (3)], the force computed above is exactly what 
those particles feel. The Sip in Eq. (23) is measured in the 
dark-matter rest frame while that in Eq. (16) is measured 



in the fundamental observer's frame. Then, in order to 
use the Sip obtained from Eq. (16) in Eq. (23), we need to 
perform a gauge transformation V ' a ip — > \7 a <p—pv a where 
v a is the peculiar velocity of the dark-matter particle rel- 
ative to the fundamental observer and consequently we 
will have an additional term -jj-pv a in Eq. (23) which is 
the velocity-dependent term identified in [27] . From here 
on we will always use the V ' a ip measured in the funda- 
mental observer's frame (in which the density perturba- 
tion is obtained more directly), and so Eq. (23) should 
be changed to 



1 bc X X 



We stress that here 



C(ip) 



(v a 



ip — ipv 



v = f = ax + ax = ax 



(24) 



(25) 



C. Code Units 

In our numerical simulation we use a modified version 
of MLAPM ([39], see III A for a brief description of the code 
and the essential modifications to it), and we will have to 
change or add our Eqs. (16, 22, 28, 29, 30). For this, the 
first step is to convert the quantities to the code units of 
MLAPM. Here, we briefly summarize the main features. 

MLAPM code uses the following internal units (where a 
subscript c stands for " code" ) : 



X c 


= x/B, 


Pc 


= P/(H B) 


tc 


= tH 


i>, 


= ^>/(H B) 2 


Pc 


= p/p, 


a 


= ac 2 y/K,Sip/ (HqB 



(31) 



where r ? is the velocity of a dark-matter particle relative 
to the fundamental observer which is at the same position 
at that moment so that x l = but x % 7^ 0. In the non- 
relativistic limit, the spatial components of Eq. (24) can 
be written as 



di 2 



= -V r J>- 



C v (ip) , 
CQp) 



rip- 



CJip) d(r/a) . 
-a — - — ip 



C(ip) 



dt 



(26) 



where t is the physical time coordinate. If we use the 
comoving coordinate x instead, then this becomes 



2-x = 

a 



1 



:V X $- 



cm 






ip + ipx)(27) 



where we have used Eq. (20). 

The canonical momentum conjugate to x is p = a 2 x 
so from the equation above we have 



dp 



dx p 

~dl ;: ^' 

CDM 



(28) 



1. 



cm 



dt 



a x C(ip) 

= — y $ 

a x C(ip) 

dt a 



V x (^ + a ip-k. 



e(lf) (V x ^ + ^p CD M 



(30) 



where B denotes the comoving size of the simulation box, 
Hq is the present Hubble constant, and p, with subscript, 
could represent the density of either CDM (/O c ,cdm) or 
baryons (/O c .b)- In the last line the quantity u is the scalar 
field perturbation Sip expressed in terms of code units and 
is new to the MLAPM code. 

In terms of u, as well as the (dimensionless) back- 
ground value of the scalar field, \fnip, some relevant quan- 
tities are expressed written in full as 



V(ip) = 



A 4 



(yfnip + 



B*my 



C(ip) 



r- 



cxp 



7 ( y/Kip + 



B 2 H% 



v^A 4 



\fap+^Lu 



l+a ' 



Ctp = 7^Kexp 



B 2 H 2 

i[V^P + 



ac 



(32) 



and the background counterparts of these quantities can 
be obtained simply by setting u — (recall that u repre- 
sents the perturbed part of the scalar field) in the above 
equations. 

Using all the above newly-defined quantities, we can 
rewrite Eqs. (16, 22, 28, 29, 30) as 



where Eq. (29) is for CDM particles and Eq. (30) is 
for baryons. Note that according to Eq. (29) the quan- 
tity a\og[C(ip)] acts as an effective potential for the fifth 
force. This is an important observation and we will re- 
turn to it later when we calculate the escape velocity of 
CDM particles within a virialized halo. 

Eqs. (16, 22, 28, 29, 30) will be used in the code to eval- 
uate the forces on the dark-matter particles and evolve 
their positions and momenta in time. 



dx c 
dt c 


Pc 

a 2 ' 




dp c 
dt c 


a 


1 c 

P^- (Vu + a^/K$p c ) 

a yJnC 


?2 $ c = 


= ^cdmC 


Y ( Pc.CDM"^ — 1 J 




3 r^ / 

-1 fin n 


^ V ~V 3 



2-— "' " H 2 



(33) 

(34) 



(35) 



and 



V 2 u 



-^n 



cdmCj; 



C v 

PcCDM-s; 1 



■Vu, - K 



V ^3, 



#o 2 



a%36) 



where Q C dm = 87rG/9cDM/3i?o an d ^b = St:Gpb/3Hq 
are the dark matter and baryonic fractional energy den- 
sities at the present time. Note that in Eq. (34) the term 
in the brackets on the right-hand side only apply to dark 
matter and not to baryons. Also note that from here on 
we shall use V = <9 Xc , V 2 = <9 Xc • 9 Xc unless otherwise 
stated, for simplicity. 
We also define 



A = 



kA 4 



(37) 



which will be used frequently below. 

Making discrete versions of the above equations for TV- 
body simulations is then straightforward, and we refer 
the interested readers to Appendix A to the whole treat- 
ment, with which we can now proceed to do TV-body sim- 
ulations. 



III. SIMULATION DETAILS 
A. The TV-Body Code 

We have modified the publicly available TV-body code 
MLAPM (Multi-Level Adaptive Particle Mesh) for our sim- 
ulations. This code uses multilevel grids [40-42] to ac- 
celerate the convergence of the (nonlinear) Gauss-Seidel 
relaxation method [41] in solving boundary value partial 
differential equations. Furthermore, it is also adaptive 
and refines the grid in regions where the mass/particle 
density exceeds a certain predefined threshold. Each re- 
finement level forms a finer grid (which might contain 
many parts that are spatially disconnected) which the 
particles will be then linked onto and on which the den- 
sities are computed, the scalar field equation and Poisson 
equation are solved, the total force on the particles are 
obtained, and the particles are drifted and kicked with a 
smaller time step. Thus MLAPM has two kinds of grids: the 
domain grid which is fixed at the beginning of a simula- 
tion, and refined grids which are generated according to 
the particle distribution and which are destroyed after a 
complete time step. 

One benefit of such a setup is that in low-density re- 
gions where the resolution requirement is not high, fewer 
time steps are needed, while the bulk of the computing 
sources can be used in the few high-density regions where 
high resolution is needed to ensure precision. 

Some technical issues must be controlled. For exam- 
ple, once a refined grid is created, the particles in that 
region will be linked to it and densities on it are cal- 
culated, then the coarse-grid values of the gravitational 
potential are interpolated to obtain the corresponding 
values on the finer grid. When the Gauss-Seidel iteration 



is performed on refined grids, the gravitational potential 
on the boundary nodes is kept constant and only those 
values on the interior nodes are updated according to 
Eq. (A6) so as to ensure consistency between coarse and 
refined grids. This point is also important in the scalar 
field simulation because, like the gravitational potential, 
the scalar field value is also evaluated on and communi- 
cated between multi-grids (note that different boundary 
conditions might lead to different solutions to the scalar 
field equation of motion) . 

In our simulation, the domain grid (the finest grid 
that is not a refined grid) has 128 3 nodes, and there 
is a ladder of coarser grids with 64 3 , 32 3 , 16 3 , 8 3 , 4 3 
nodes respectively. These grids arc used for the multi- 
grid acceleration of convergence: for the Gauss-Seidel re- 
laxation method, the convergence rate is high for the 
first several iterations, but then quickly becomes very 
slow; this is because the convergence is only efficient for 
the high-frequency (short-range) Fourier modes, while 
for low-frequency (long-range) modes more iterations do 
not help much. To accelerate the solution process, one 
then switches to the next coarser grid - for which the 
low-frequency modes of the finer grid are actually high- 
frequency ones - and the convergence is speeded up. The 
MLAPM Poisson solver adopts the self-adaptive scheme: if 
convergence becomes slow on a grid, then go to the next 
coarser grid where it is expected to be faster, if conver- 
gence is achieved on a grid, then interpolate the relevant 
quantities back to the finer grid (provided that the latter 
is not on the refinements) and solve the equation there 
again. This goes on indefinitely until a converged solution 
on the domain grid is obtained, or until one arrives at 
the coarsest grid (normally with 2 3 nodes) on which the 
equations can be solved exactly using other techniques, 
or by simply iterating many times until convergence is 
achieved. For the scalar-field equation of motion, we find 
that with the self-adaptive scheme in certain regimes the 
nonlinear Gauss-Seidel solver tends to fall into oscilla- 
tions between the coarser and finer grids; to avoid such 
situations, we then use V-cycle [41] instead. 

For the refined grids the method is different. Here, we 
just iterate Eq. (A6) until convergence, without resorting 
to coarser grids for acceleration. 

As is normal in the Gauss-Seidel relaxation method, 
convergence is deemed to be achieved when the numerical 
solution u\, after n iterations on grid k, satisfies that the 
condition that the norm || • || (mean or maximum value 
on a grid) of the residual 

e k = L k (u k n )~f k , (38) 

is smaller than the norm of the truncation error 

r k = L k ~ l {1lul)-1l[L k {u k n )\ (39) 

by a certain amount, or, in the V-cycle case, that the 
reduction of residual after a full cycle becomes smaller 
than a predefined threshold (indeed the former is satis- 
fied whenever the latter is). Note here that L k is the dis- 
cretization of the differential operator Eq. (A4) on grid k 



and L h ~ x a similar discretization on grid k — 1, fk is the 
source term, 1Z is the restriction operator to interpolate 
values from the grid k to the grid k — 1. In the modified 
code we have used the full-weighting restriction for 1Z. 
Correspondingly, there is a prolongation operator V to 
obtain values from grid k — 1 to grid k, and we use a 
bilinear interpolation for it. For more details see [39]. 

MLAPM calculates the gravitational forces on particles 
by a centered difference of the potential $ and inter- 
polates the forces at the locations of particles by the 
so-called triangular-shaped-cloud (TSC) scheme to ensure 
momentum conservation on the grid. The TSC scheme is 
also used in the density assignment given the particle 
distribution. 

Some of our main modifications to the MLAPM code for 
the coupled scalar field model are: 

1. We have added a parallel solver for the scalar field, 
based on Eq. (A2). It uses a nonlinear Gauss-Seidel 
scheme for the relaxation iteration and the same 
criterion for convergence as the default Poisson 
solver. But is adopts a V-cyclc instead of the self- 
adaptive scheme in arranging the Gauss-Seidel it- 
erations. 

2. The value of u solved thereby is then used to calcu- 
late the total matter density, which completes the 
calculation of the source term for the Poisson equa- 
tion. The latter is then solved using fast Fourier 
transform on the domain grids and self-adaptive 
Gauss-Seidel iteration on refinements. 

3. The fifth force is obtained by differentiating the u 
as the gravity is computed by differentiating the 
gravitational potential. 

4. The momenta and positions of particles are then 
updated, or in other words the particles are kicked 
and drifted, where in the kicks we take into account 
both gravity and the fifth force. 

There are a lot of additions and modifications to ensure 
smooth interface and the newly added data structures. 
For the output, as there are multilevel grids all of which 
host particles, the composite grid is inhomogeneous and 
so we choose to output the positions and momenta of the 
particles, plus the gravity, fifth force and scalar field val- 
ues at the positions of these particles. We can of course 
easily read these data into the code, calculate the cor- 
responding quantities on each grid and output them if 
needed. We also output the potential and scalar field val- 
ues on the 128 3 domain grid. 



B. Physical and Simulation Parameters 

The physical parameters we use in the simulations are 
as follows: the present-day dark-energy fractional energy 
density f^DE = 0.743 and Q m = Ocdm + ^b = 0.257, 



H a = 71.9 km/s/Mpc, n s = 0.963, er 8 = 0.761. Our 
simulation box has a size of 64/i _1 Mpc, where h = 
.ffo/(100 km/s/Mpc). We simulate 4 models, with pa- 
rameters a = 0.1 and 7 = -0.05,-0.10,-0.15,-0.20 
respectively (such parameters are chosen so that the de- 
viation from ACDM will not become too large to be re- 
alistic). In all those simulations, the mass resolution is 
1.114 x 10 9 /i _1 Mq; the particle (both dark matter and 
baryons) number is 256 3 (see Sect. Ill C for a discussion); 
the domain grid is a 128 x 128 x 128 cubic and the finest 
refined grids have 16384 cells on each side, corresponding 
to a force resolution of about 12/i _1 kpc. 

We also make a run for the ACDM model using the 
same physical parameters and different initial condition 
(see Sect. HIE). 



C. Distinguishing Baryons from Dark Matter 

In our simulations only dark matter particles are cou- 
pled to the scalar field. The baryons do not contribute 
to the scalar field equation of motion, nor are they in- 
fluenced by the scalar-field fifth force. Therefore, it is 
important to make sure that they are distinguished and 
treated appropriately. 

In the modified code we distinguish baryons and dark 
matter particles by tagging them differently. We consider 
the situation where 17.12% of all our matter particles are 
baryonic (fi B = 0.044) and 82.88% arc dark matter l . 
During the generation of the initial condition (initial dis- 
tribution and displacements of particles), we loop over 
all particles and for each particle we generate a random 
number from a uniform distribution in [0, 1]. If this ran- 
dom number is less than 0.1712 then we tag the particle 
as baryon, otherwise we tag it as dark matter. Once these 
tags have been set up they are never changed again, and 
the code then determines whether or not a particle con- 
tributes to the scalar field evolution and feels the fifth 
force according to its tag. 

Often in A^-body simulations, the number of particles 
is chosen to be the same as the number of cells in the 
domain grid, but in our simulations we have set the 
former (256 3 ) to be eight times bigger than the latter 
(128 3 ). This choice obviously increases the mass resolu- 
tion, but what is more important for us is that it produces 
a smoother and more reasonable dark matter density on 
the grid. In order to see this, remember that only 82.88% 
of all particles are dark matter, which means that, if we 



1 It is time to stress that our Qqom = SttGpcdm/^Hq is the 
fractional energy density of the bare dark matter particles, which 
is not weighed by the coupling function C(ip). As we mentioned 
in Appendix C and [17], Pcdm c< a~ 3 as in ACDM but the 
behaviour of C((fi)pcDM (which is the actual quantity appearing 
in the Poisson equation) can be rather complicated. In all models 
we are simulating, including the ACDM one, we have the same 
PCDMi rather than C(</j)pcDM> a t present-day. 
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FIG. 1: (Color Online) Figures to illustrate the background evolution of our coupled scalar field models. Upper Left Panel: The 
fractional energy densities for radiation (green), coupled dark matter (black), baryons (blue) and the scalar field (red); note that 
ficc = 8ivGC (ip) pcdm /SH 2 . Upper Right Panel: The equation of state of the scalar field, w = p v /p v in which p v = ^ip 2 — V(ip) 
and p v — 7;ip 2 + V(ip). Lower Left Panel: The ratio between the Hubble expansion rate in the coupled scalar field model and 
that in the ACDM paradigm, other physical parameters such as Pb, Prad, Pcdm being held the same, as a function of the scale 
factor a. Lower Right Panel: The "varying mass" of the dark matter particles as a function of a - here mo is the constant bare 
mass of the particles and m = e 7V mo- In all figures we have chosen a = 0.1 and the solid, dotted, dashed and dot-dashed curves 
represent the models with 7 = —0.05, —0.10, —0.15 and —0.20 respectively. 



use 128 3 rather than 256 3 particles, about 1/5 of all the 
grid cells do not contain dark matter particles at the 
initial time. The resulting dark matter density therefore 
shows some artificial discreteness, which not only does 
not reflect reality but might also cause the solver for the 
scalar field equation to diverge. By having more parti- 
cles, we can make the dark matter densities smoother 
and resolve the problem of over-discreteness. 



D. Background and Linear Perturbation Evolution 

In general, a coupling between the scalar field and the 
dark matter particles not only affects the force law and 
clustering properties of those particles [as described by 



Eqs. (16, 22, 28, 29, 30)], but it also influences the back- 
ground and linear perturbation evolution of the Universe. 
The modification in the background cosmic expansion 
rate directly changes the rate of the matter clustering, 
while the modification in the linear perturbation growth 
might lead to a different initial condition for the N- 
body simulation from the ACDM result. Consequently, 
we must take appropriate care of these issues in the N- 
body simulations in order to obtain reliable and complete 
numerical results. 

The background expansion rate of the universe is com- 
pletely governed by (the zeroth order parts of) the scalar 
field equation of motion 



(p + 3Hip + V V + C^pcdm 
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FIG. 2: (Color Online) The linear power spectra of the coupled scalar field model. Upper Left Panel: The Cosmic Microwave 
Background (CMB) power spectra for the four models with a — 0.1. Upper Right Panel: The same but for the models with 
a = 0.5. Lower Left Panel: The matter power spectra at present day (redshift 2 = 0) for the four models with a = 0.1. Lower 
Right Panel: The same but for the models with a = 0.5. In all figures the solid (black), dotted (blue), dashed (green) and 
dot-dashed (red) curves represent the models with 7 = —0.05, —0.10, —0.15 and —0.20 respectively. 



the Friedman equation 
3H 2 = k 



PB + PRAD + C{lfi)pcDM + X^ 2 + V(<f) |(41) 



where H = a/a is the Hubble expansion rate and prad 
includes contributions from photons and masslcss neutri- 
nos (i.e., radiation), and the Raychaudhuri equation 



3(. H 

K r 

~2 



H* 

PB + 2/9RAD + C(<f)p C DM + 2(p 2 ~ 2V(<p)] (42) 



The introduction of Eq. (40) introduces a new degree 
of freedom ip to the system, for which the initial condi- 
tion and relevant parameter (the A in V(ip)) should be 
chosen appropriately to guarantee consistency. More ex- 
plicitly, we have to adjust the values of A and '/'im , ^ini 



so that today H = Hq where Hq is the measured Hubble 
constant which we set to be 71.9 km/s/Mpc in this work. 
This is a nontrivial requirement and is achieved through 
a trial-and-error process. We shall leave all the details of 
the numerical algorithm to Appendix C, as they are not 
our major concern here. 

In Fig. 1 we have shown some representative results 
for the influences of the scalar field coupling on the back- 
ground cosmology. Clearly, for fixed a, increasing \j\ 
means increasing the coupling strength, leading to more 
dramatic evolution of the scalar field, and this is why in 
the lower right panel we can see that as |— y | increases, so 
the quantity e lv falls increasingly below 1. Meanwhile, 
as the scalar field rolls faster for larger \y\, the equation 
of state parameter w approaches —1 (potential-energy- 
dominated regime) later in time (the upper right panel). 
Also, because e lv deviates more from 1 , the contribution 
of the coupled dark matter (p cc = /ccdmc 7 ^ in which 
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FIG. 3: (Color Online) Upper Panels: The matter power spectra for baryons only {Left) and for dark matter only {Right), 
at current time (2 = 0). Lower Panels: The same but at an earlier time z — 49 where the initial condition for the A^-body 
simulations are computed. In all figures a = 0.1 and the solid (black), dotted (blue), dashed (green) and dot-dashed (red) 
curves represent the models with 7 = —0.05, —0.10, —0.15 and —0.20 respectively. 



Pcdm oc a~ 3 , see Appendix C) to the Fricdmann equa- 
tion [Eq. (41)] is smaller and the contributions from other 
matter species are greater (upper left panel). Because 
dark matter is the principal ingredient driving the ex- 
pansion of the Universe in the matter-dominated era, this 
in turn implies that the expansion rate will differ from 
the ACDM prediction (lower left panel). Note that the 
change of expansion rate at early times necessarily cause 
the matter power spectrum at z ~ 50 to differ from the 
ACDM result, a point we will return to shortly. 



low i. The matter power spectrum, on the other hand, 
is increased on all scales: on large scales this is purely 
because the universe is now expanding more slowly, al- 
lowing matter to cluster more. On small scales there is 
another effect - the boost due to the scalar field fifth 
force, which is almost always proportional to gravity in 
magnitude and parallel in direction (see below). It is clear 
from the figure that the results are rather insensitive to 
the parameter a, and so in what follows we shall only 
consider the case of a = 0.1. 



In Fig. 2 we have displayed the effects of the scalar 
field coupling on the linear cosmic microwave background 
(CMB) and matter power spectra. Because of the de- 
crease in the expansion rate, the angular diameter dis- 
tance increases so that the peaks of the CMB spectrum 
are shifted rightwards, and on large scales the pertur- 
bation in the scalar field (cf. Fig. 4 below) modifies the 
integrated Sachs- Wolfe effect to increase the power at 



Since wc need to know the matter power spectrum at 
early time z ~ 50 in order to general initial conditions 
for the A-body simulations, we also plot it in Fig. 3. We 
further separate the baryons from dark matter, the latter 
being our major concern: on small scales it is clear that 
the dark matter power spectrum is greater than that of 
baryons, due to the extra fifth force it experiences (cf. up- 
per panels). Most interestingly, we see that even at z — 49 
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FIG. 4: Upper Left Panel: The fractional change of the growth factor D+ of the dark matter density perturbation in the coupled 
scalar field model as compared with the ACDM prediction; for clearness shown are only the results for the two models with 
a = 0.1,7 = —0.05 and a — 0.1,7 = —0.15, as indicated above the curves, and for each model the solid, dotted, dashed and 
dot-dashed curves represent respectively the result for k = 0.0001, 0.001, 0.01 and 0.1 Mpc -1 . Upper Right Panel: The same as 
above, but for the fractional change of D+ as a function of the scalar factor a. Lower Left Panel: The time evolution of the 
fractional energy density of the kinetic energy of the scalar field, fikin = ^ivGip 2 /3H 2 for the two models with a = 0.1, 7 = 0.05 
(solid curve) and a — 0.1,7 = —0.15 (dotted curve) respectively. Lower Right Panel: The evolution of the density contrast of 
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the dark matter power spectra for different models can be 
rather different, again due to the modified expansion rate 
and the fifth force. This means that in our iV-body sim- 
ulations we should not use the same initial condition (as 
in [17, 19], where the chameleon effect is so strong that at 
early times the matter power spectrum is indistinguish- 
able from that of ACDM). Instead, we should generate 
initial conditions separately for different models (in our 
case different values for 7) - this is the topic of Sect. Ill E 
and Appendix A. 



E. Initial Conditions for the A-Body Code 

As we described in detail above, the fact that the 
scalar-field-dark-matter coupling begins to take effect at 
rather high redshift indicates that the initial conditions 
for the N-body simulations of our coupled scalar field 
models (at redshift zi ~ 50 here) will also be different 
from those in ACDM. To account for this modification, 
the most straightforward approach is to generate the lin- 
ear matter power spectra for our coupled scalar field mod- 
els at the starting redshift Zj ~ 50, and utilize them to 
produce the Gaussian random density fluctuation field 
and displace the particles. In GRAFIC2 the matter power 
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FIG. 5: (Color Online) Snapshots of the particle distribution in our four coupled scalar field models as indicated by the subtitles 
of the panels, a is the scale factor and a = 1 is the present time. For clearness we only pick out a slice of the simulation box 
with 30 /i _1 Mpc < z < 30.3 /i -1 Mpc and /i _1 Mpc < x,y < 64 /i -1 Mpc. The blue dots represent dark matter particles and 
red dots baryons. 



spectrum is generated at time z = 0, normalized to a pre- 
selected value and then evolved back to Zi and used to 
displace particles. We shall not follow this in our work; in- 
stead we adopt the same initial condition in CAMB at red- 
shift z ~ 10 6 for all models including the ACDM one and 
evaluate the power spectrum at zi - in this way the linear 
08 today will be different for different models (in contrast 
erg is the same for all models in some other works). This 
is because we are interested in how the scalar field affects 
the evolution of structure as compared to ACDM, given 
the same initial condition at very early times. 

In principle, the matter power spectrum at z% that 
is used in GRAFIC2 needs to be generated for each 
model, for example by linear perturbation code that 
is default in GRAFIC2; however, here wc adopt a 
more economical method. To sec this, in the up- 
per panels of Fig. 4 we have displayed the time 
evolution of [D + — -D+,acdm] /-D+.acdm as well as 

D+ - -D+,acdm /-D+,acdm, with D + being the linear 
growth factor we have discussed in Appendix D. For 
simplicity, we only show these for the two models with 
a = 0.1 and 7 = —0.05,-0.15 as indicated in the fig- 
ure, and for each model the solid, dotted, dashed and 
dash-dotted curves represent respectively the results for 



k = 0.0001, 0.001, 0.01 and O.lMpc . There are several 
important features in these plots. First, at early times 
(z > 10 4 ) we see that the difference is negligible as the 
scalar field has yet to take effect. Second, on large scales 
(small k) the results tend to converge to a curve which 
describes how the change of background expansion rate 
modifies the density growth (the fifth force has no ef- 
fect here because the scale is out of its range). Third, on 
smaller scales (large k) the fifth force begins to act and 
further enhances the growth of dark matter density per- 
turbation. These plots show that D+, the dark matter 
linear growth factor, depends on both k and time as we 
mentioned in Appendix D, and this must be taken into 
account. 

As detailed in Appendix D, the displacement and pe- 
culiar velocity of particles are given by D + d and D + d 
respectively, where d is a vector which is the same for 
different models (remember again wc start with the same 
initial condition in CAMB at z ~ 10 6 ). Consequently, the 
differences in the displacements and peculiar velocities 
in our models are equivalent to the differences in D + (zi) 
and -D+ (z{). Instead of evaluating D + {zi) and D + (zi) 
with (a modified version of) GRAFIC2 for each model, we 
use the ACDM results for D + (zj) and D + (zi) computed 
by GRAFIC2, and also calculate /j,\ = Z? + /D +] acdm and 
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FIG. 6: (Color Online) The relation between the magnitudes of the fifth force and that of gravity for the four coupled scalar field 
models at three different output times a — 0.2, 0.5, 1.0. The black solid line in each panel represents the analytical approximation 
/ = 2j 2 ■ ^ DM g (see text) and the ~ 150, 000 green dots the results from the simulations. The particles are the same as those 
in Fig. 5. 



fj,2 = -D+/-D+,acdm at Zi using CAMB. Then the particle 
displacement and peculiar velocity in the coupled scalar 
field models are simply given respectively by \i\^ times 
those generated by GRAFIC2. 

Furthermore, we shall take into account the differences 
between the coupling and non-coupling matter species as 
following. The baryons do not feel the fifth force (remem- 
ber that 'baryons' here simply means non-coupling dark 
matter particles, and this is why we do not use the baryon 
matter power spectrum generated by CAMB to displace 
these particles). The difference between the baryonic lin- 
ear growth factor D + from -D+,acdm is mainly due to the 
modified expansion rate, and so we use the D + calculated 
for very small k to evaluate \i\^. for baryons. Dark matter 
particles, on the other hand, do feel the fifth force which, 
on small scales, has a magnitude of l^ 2 times that of 
gravity, and so we use the D+ calculated for large k to 
evaluate /ii.2 for dark matter. In this way, the bias be- 
tween the two matter species is approximately reflected 
in our initial conditions. The scale dependence of D + 
does not appear to have been considered in generating 
initial conditions in the earlier literature. 

When generating the initial displacement and peculiar 
velocities of particles, we run a random number genera- 
tor for a variable with a uniform distribution U[0, 1]. As 



mentioned above, if the number generated is less than 
^b/ (^b + ^cdm), we tag the particle as a baryon and 
use /j,i t 2 to displace it; otherwise we tag the particle as 
dark matter and use jx\^ to compute its displacement 
and velocity. Once this has been done, no particle will 
ever be re-tagged so that the consistency is not spoiled. 
Those particles tagged as dark matter will contribute to 
the scalar field evolution and be influenced by the scalar 
field fifth force, while those which are tagged as baryons 
will not. 



IV. SIMULATION RESULTS 

In this section we present the main results from the N- 
body simulations we have performed. We start by listing 
some preliminary results to give a rough idea about the 
effects of the couplings we are studying. 



A. Preliminary Results 

In Fig. 5 we show snapshots of the particle distribution 
at different output times. We can see clearly the trend of 
matter clustering. The blue and red dots are the dark 
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FIG. 7: (Color Online) The gravitational potential on the z plane where z = 32/i _1 Mpc. Dark regions are where the potential 
is deeper while light regions are where it is shallower. The four columns are for the four models we consider and the three rows 
are for three output times, respectively a = 0.2, 0.5 and 1 where a is the cosmic scale factor. 



matter and baryons particles respectively. 

Since one of the most important influences of the scalar 
field coupling is to exert a fifth force on the dark matter 
particles, we arc interested in the magnitude of the fifth 
force as compared with that of gravity. This is given in 
Fig. 6, where we have again displayed the results for the 
simulated models at different output times. To under- 
stand this figure, note that if the contributions from the 
scalar field potential to the Poisson equation and scalar 
field equation of motion are negligible (which turns out 
to be the case in our simulations) and if all particles 
(dark matter and baryons) have the same coupling to 
the scalar field, then the right hand sides of Eqs. (35, 36) 
are proportional, with a coefficient 2j, which implies that 
u = 27$ c ; thus Eq. (34) says that the strengths of fifth 
force (/) and gravity (g) should satisfy / = 2-f 2 g. 

In our models only a fraction $7cDM/^m> where Cl m = 
f^CDM + ^Bi of all the particles are coupled to the scalar 
field. Thus, assuming dark matter and baryonic particles 
arc distributed in the same way, we should have 



/ = 2 7 



2^CDM 



- -/ 



-9 



(43) 



because only a fraction f^cDM/^m of the particles which 
produce gravity also produce the fifth force. Eq. (43) 
is plotted in Fig. (i as the solid lines. The dots are the 



simulation results for f/g for the particles outputted in 
Fig. 5. The agreement between the analytical approxima- 
tion and numerical solution is remarkably good, which 
serves as a test of our Newton-Gauss-Seidel solver of the 
scalar field equation of motion 2 . 

Furthermore, Fig. 6 also shows some clearly different 
feature from the results of [19], in which the fifth force is 
suppressed by the chameleon mechanism in certain cases 
(such as at early times and in high-density regions). Here, 
the fifth force is not suppressed because the potential 
V(<p) is negligible, and the approximation Eq. (43) works 
well anywhere all the time. This confirms our earlier claim 
that the matter power spectrum will be changed at early 
times as will be the initial condition for the iV-body code. 

As should evident now, the spatial configuration of the 
scalar field ip, or equivalently the rescaled scalar field u, 
should closely follow that of the gravitational potential 
$ c , for the same reasons that led to Eq. (43). Figs. 7 and 
8 confirm this: as can be seen there, for all the models 
and all the output times, the configurations of — u 3 and 



2 Remember that the Poisson equation on cither the domain grid 
or the refinement is solved using different methods from those 
for the scalar field equation of motion. 
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FIG. 8: (Color Online) The same as Fig. 7 but for — ip where <p is the value of the scalar field. The negative sign is added to 
make the plot look similar to Fig. 7, and otherwise they will be compensating for each other. 



(f> c arc indistinguishable. 

The above figures show that in our coupled scalar field 
models the fifth force / (scalar field u) closely mimics the 
gravity g (gravitational potential $ c ), and so the approxi- 
mation / = 2^/ 2 g, used in many other previous simulation 
works {e.g. [27]), is a good one. 



B. Matter Power Spectrum 

The nonlinear matter power spectrum measured from 
the output particle distribution is of more observational 
interest and it is the theme of this subsection. The mat- 
ter power spectrum in the present work is measured us- 
ing POWMES [43], which is a publicly available code based 
on the Taylor expansion of trigonometric functions and 
yields Fourier modes from a number of fast Fourier trans- 
forms controlled by the order of the expansion. 

Fig. 9 displays the fractional changes of the matter 
power spectra P(k) due to the scalar field coupling, where 
AP = Pgcalar — -Pacdm ■ For comparison, we have also 



3 Note we expect that u = 27<E> C approximately as explained above, 
and 7 < 0: the minus sign here ensures that Fig. 8 follows Fig. 7, 
rather than compensating it, so that we have a better visual 
impression. 



shown the linear results as dashed curves. We can see that 
on large scales (small k) there is an overall agreement be- 
tween the linear and nonlinear predictions, which is not 
surprising since those scales have not experiences much 
nonlinear complication. The agreement is very good at 
early times (the top panels) when nonlinear effects gen- 
erally have not taken place on scales k ~ 0.1 /iMpc ; at 
late times (top curves in bottom panels) there is a slight 
mismatch, which is understandable because our simula- 
tion box is not big enough to allow measurement of P(k) 
below k = 0.1 /iMpc^ 1 where nonlinear effects already 
enter. Indeed, even in the latter case, we can see the 
trend that the linear and nonlinear curves merge towards 
at small-fc. The fact that the nonlinear P(k) reduces to 
the linear one on large scales is not trivial as the back- 
ground cosmology in our TV-body simulations has also 
been modified 4 (as it is in CAMB, which is used for the 
linear calculation), and we take this as another test of 
our modified MLAPM code. 

The main advantage of the TV-body simulation is its 
ability to probe the nonlinear structure formation and so 



4 As mentioned above: the larger | — y | , the slower the universe ex- 
pands (cf. Fig. 1, lower left panel), and thus the more growth the 
matter perturbations experience and the larger P(k) becomes. 
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FIG. 9: (Color Online) The linear and nonlinear matter power spectra P(k) as compared with those of the ACDM paradigm. 
Shown are the fractional change of P(k) for the four models with a = 0.1 and 7 = —0.05 (black curves), —0.10 (green), —0.15 
(blue), —0.20 (purple) and at four different output times a = 0.3 (Upper Left Panel), 0.5 (Upper Right), 0.7 (Lower Left), 1.0 
(Lower Right). Solid curves are from iV-body simulations and dashed curves from linear perturbation calculation. 



we are more concerned with the P(k) on smaller scales, 
which is also plotted in Fig. 9. We can see that on in- 
termediate scales (k ~ 1 /iMpc -1 ) the nonlinear P(k) 
beats the linear one, but on even smaller scales it falls 
behind the linear value again. For the four models in our 
simulations the enhancements of the matter power range 
from negligible to ~ 100%, but much of the enhancement 
is due to the modified background expansion and is al- 
ready seen in the linear results. More interestingly, the 
enhancement is more significant at earlier times. This can 
been understood as follows: at later times the contribu- 
tion of the dark matter density to the Poisson equation 
(and thus the gravitational potential) is decreased due 
to the coupling function C(ip) becoming increasingly less 
than 1, weakening further clustering of the (not only dark 
but also baryonic) matter. Note that such a trend can 



also be seen in the linear P(k), by comparing the bottom 
panels, although it is not so obvious. 

In order to display the bias developed between baryons 
and dark matter, we have also plotted the P(k) for 
baryons and dark matter separately in Fig. 10. As ex- 
pected, the P(k) for dark matter is always larger, signi- 
fying a stronger clustering due to the assistance from the 
fifth force. The larger I7I is, the stronger the fifth force 
will be, and the larger the bias will become. 



C. Mass Function 

We identify halos in our iV-body simulations using MHF 
(MLAPM Halo Finder) [44], which is the default halo finder 
for MLAPM. MHF optimally utilizes the refinement structure 
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FIG. 10: (Color Online) The bias between the nonlinear matter power spectra P(k) of dark matter and baryons for the four 
models we are simulating, as indicated in the subtitles of the panels. The solid curves are for the output time a = 0.5 and 
dashed curves for a — 1.0. For each style of curve the lower (black) one is for baryons and the upper (purple) one is for dark 
matter. 



of the simulation grids to pin down the regions where po- 
tential halos reside and organize the refinement hierarchy 
into a tree structure. MLAPM refines grids according to the 
particle density on them and so the boundaries of the 
refinements are simply isodensity contours. MHF collects 
the particles within these isodensity contours (as well as 
some particles outside). It then performs the following 
operations: (i) assuming spherical symmetry of the halo, 
calculate the escape velocity v esc at the position of each 
particle, (ii) if the velocity of the particle exceeds v esc 
then it does not belong to the virialized halo and is re- 
moved. Steps (i) and (ii) are then iterated until all un- 
bound particles are removed from the halo or the number 
of particles in the halo falls below a pre-defined threshold, 
which is 20 in our simulations. Note that the removal of 
unbound particles is not used in some halo finders using 



the spherical ovcrdcnsity (SO) algorithm, which includes 
the particles in the halo as long as they are within the 
radius of a virial density contrast. Another advantage of 
MHF is that it does not require a pre-defined linking length 
in finding halos, such as the friend-of-friend procedure. 

As explained in detail in [19], part of the MHF algo- 
rithm also needs to be modified for the coupled scalar 
field model, because the scalar field ip behaves as an ex- 
tra "potential" (which produces the fifth force) and so the 
dark matter particles experience a deeper total "gravita- 
tional" potential than in the ACDM paradigm, all other 
things being equal. Consequently, the escape velocity for 
dark matter particles increases compared with the New- 
tonian prediction. As the dark matter particles in the 
coupled scalar field simulations are typically faster than 
in the ACDM simulation, if we underestimate v esc then 
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some particles which should have remained in the halo 
will be incorrectly removed by MHF. In [19] it was found 
that the underestimate of the mass function by using de- 
fault MHF in the coupled scalar field models could be up 
to a few percent. 

In the models we are considering here, things are even 
more complicated. In [19], the chameleon mechanism en- 
sures that the scalar field takes a very tiny value every- 
where, all the time (typically y/Hip < 10 -5 ) and so the 



coupling function C((p) = e 



■Jy/Klfl 



1 is a very good ap- 



proximation - this means that the contribution of the 
dark matter to the Poisson equation is not significantly 
modified compared with that in ACDM. In the models 
here, however, C(<p) can be up to 30% less than 1, so that 
the source term of the Poisson equation is significantly 
different from ACDM. 

Obviously, we need to take both of these two major 
changes into account when we design a new algorithm 
to identify halos. An exact analytical calculation of the 
escape velocity and Poisson source term in the coupled 
scalar field model turns out be difficult, and so we intro- 
duce an approximate alternative, which is based on the 
MHF default method [45]. 

The default MHF code works out v esc using the Newto- 
nian result 



2|$l 



(44) 



in which $ is the gravitational potential. Under the as- 
sumption of spherical symmetry for the halos, the Poisson 
equation V 2 $ = 4nGp m could be integrated once to give 



dr 



GM(< r) 



(45) 



which is just the Newtonian force law. This equation can 
be integrated once again to obtain 



$(r) 



G 



M(< r') 



dr' + $ 



(46) 



where $o is an integration constant and can be fixed [45] 
by requiring that $(r = oo) = as 



<I>, 



GM V 



G 



M(< r') 



dr' 



(47) 



in which R v i r is the virial radius of the halo and M V i r is 
the mass enclosed in R V i r . 

In the coupled scalar field models here, the two modi- 
fications mentioned above are reflected in Eqs. (44, 45). 

For Eq. (44), we realize that the gravitational potential 
$ is not the only factor determining the escape velocity: 
there is also a contribution from the scalar field (p. There 
are different ways to approach this problem, the simplest 
of which is to obtain the total "potential" by rescaling $ 
with a constant (remember we have shown in Sect. IV A 
that ip ex $ is a good approximation all the time and ev- 
erywhere). A more delicate approach can be devised by 
noting that in the default MHF code, Eq. (45) is used in 



the numerical integrations to obtain both $(r) and $o 
[cf. Eqs. (46, 47)]. More explicitly, the code loops over all 
particles in the halo in ascending order of the distance 
from the halo centre, and whenever a particle is encoun- 
tered its mass is uniformly distributed into the spherical 
shell between the particle and its previous particle (the 
thickness of the shell is now the dr of the integration). 
When the fifth force is added, we call its contribution to 
the total "gravitational" potential $ v and its value at 
infinity, $^0- 0m" estimate of the escape velocity is now, 
instead of Eq. (44), given by 



2|$ + <Iv-$o-$ 



<p0\ 



(48) 



We have recorded the components of gravity and the 
fifth force for each particle in the simulation, and so the 
ratio f / g can be computed at the position of each parti- 
cle, which gives ($ v — & v o) / ($ — $o) a t that position, 
which can be used in Eq. (48). In this way we have, at 
least approximately, taken into account the fifth force 
and <J>ip. 

For Eq. (45), which is used to compute <f> and $o, we 
know that it is applicable in the coupled scalar field model 
as well, because the contribution from the scalar field 
potential is negligible. But one should be careful when 
interpreting the mass M here: it is not the total bare 
mass of all the particles (Mb + Mcdm) within radius 
r, but Mb + C(</j)Mcdm- Such a fact can be taken into 
account by multiplying the mass of each dark matter (not 
baryonic) particle by C((p), which is evaluated using the 
ip at the position of that particle (which we have recorded 
too). Note that observations of the masses are made by 
registering gravitational effect, and so the halo mass we 
measure should be Mb + C(<p)Mcdm rather than Mb + 
Mcdm^ in what follows we will always talk about the 
former unless otherwise stated. 

In Fig. 11 we have plotted the mass functions of our 
four coupled scalar field models at the current time, as 
compared with the prediction of the ACDM paradigm. 
Obviously, as the coupling between dark matter and the 
scalar field becomes stronger, more heavy halos will be 
produced during structure formation. However, it is not 
yet clear which physical effect is mainly responsible for 
such a pattern, and several factors could be crucial. For 
examples, as I7I increases, the fifth force strengthens and 
the background expansion gets slower, both of which 
would help boost the clustering of matter; on the other 
hand, the dark matter's contribution to the Poisson equa- 
tion is weakened due to C{(p) becoming increasingly less 
than 1 - this will weaken the clustering of matter. De- 
tailed study of the significance of these different physical 
effects is beyond the scope of this work. 



D. Halo Properties 

Finally, we are also interested in the properties of the 
dark matter halos in the coupled scalar field models. In 
the case of a chameleon-like scalar field, rcf [19] shows 
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FIG. 11: The mass functions at a — 1 for the four coupled scalar field models under investigation, as indicated by the subtitles 
of the panels. The dashed curves are the corresponding result of ACDM. 



that there is a strong environmental dependence on den- 
sity profile and baryon-dark-matter bias inside the halo, 
which is controlled by the model parameters in a compli- 
cated way. In the models generated here, the fifth force is 
everywhere unsuppressed and we do not expect such an 
environmental dependence as seen in [19]. However, we 
have several factors mentioned in Sect. IV C, which can 
potentially lead to interesting new features. 

We will look first at the internal density profiles, which 
can be measured down to small radii thanks to the higher 
spatial resolution of the self-refinement code. We have se- 
lected two typical halos from each simulation. Halo I is 
centred on (x,y,z) ~ (32.4, 31.5, 61. 2)h~ l Mpc, which 
is slightly different for different simulations, and has a 
virial mass M V i r ~ 1.29 x 10 14 /i _1 Mq; halo II is cen- 
tred on (x, y, z) ~ (54.8, 39.8, 35.7)/i _1 Mpc, which is also 
slightly different for different simulations, and has a virial 



M„ 



1.90 x 10 13 /i _1 Mq (note here that the 



virial masses are for the ACDM simulations; for scalar- 
field simulations they can be, and generally are, larger). 

Fig. 1 2 summarizes what we have, and it shows that 
the density profile for the coupled scalar field model is 
rather similar to that of ACDM, except that for the inner- 
most part the overdensity for the scalar field model (solid 
curves) drops below that of ACDM (dashed curves). 
Such a phenomenon has previously been detected by 
the authors of [28] and explained as mainly due to the 
velocity-dependent acceleration of dark matter particles 
[the last term in Eq. (26)]. However, here we find here 
that such a trend does not continue with increasing \j\ 
and is indeed reversed for large I7I values. To see this 
more explicitly, in Fig. 13 we have shown the quantity 
Ap(< -R)/pacdm(< R) (where p(< R) is the average 
density inside radius R and A means the difference be- 
tween the coupled scalar field model and ACDM) aver- 
aged over 10 of the heaviest halos of our simulations. This 
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FIG. 12: (Color Online) The internal density profiled of two selected halos (see text) at a = 1. The four panels are for the four 
models we have simulated, as illustrated by the subtitles. In all panels solid curves are for the coupled scalar field model and 
dashed curves for ACDM. The upper (black) curves are result for halo I and lower (green) curves for halo II. 



shows clearly that for small | — ^ | values the inner density 
is systematically lower in our coupled scalar field mod- 
els than in ACDM while this is not necessarily true for 
the density in the outer part of the halos. Also, as I7I 
increases, the lowering of inner density becomes less sig- 
nificant, and finally for 7 = —0.20 we see the opposite 
trend emerge. The physical explanations of the observed 
patterns are not of prime interest to our study here. 

Next, we can look at the bias between baryons and 
dark matter in the halos, which is displayed in Fig. 14. 
As expected, we find that the dark matter density is con- 
stantly higher than that of baryons, due to the extra force 
they feel (note that other factors, such as the cosmic ex- 
pansion rate and the modification of the source term in 
the Poisson equation, have the same influence on dark 
matter and baryons). In general, the bias increases with 
the coupling strength \~/\ which is easy to understand. 



V. SUMMARY AND DISCUSSION 

Couplings between any scalar field and (some of) the 
matter species can affect cosmology in various ways. 
There are two principal effects: modified source terms in 
the gravitational field equations (channel I) and direct 
new interactions between particles of the coupled matter 
species (channel II). 

The effects through channel II arise in an inhomo- 
gencous universe, but not in background cosmology. In 
principle, there is only effect of this sort - the fifth force, 
due to the exchanges of scalar quanta between matter 
particles. In practice, (e.g., in TV-body simulations), the 
scalar field tp is computed in the fundamental observer's 
frame while the fifth force must be evaluated in indi- 
vidual particle frames; consequently, the fifth force is 
split into two parts: the part from the spatial gradient of 
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FIG. 13: The fractional change of the internal density profile due to the dark matter being coupled to the scalar field, [p(< 
R) — /9acdm(< -R)]/pacdm(< R)- The solid curves are the average of ten of the heaviest halos in our simulation boxes and the 
dashed curve is 0. We show this for the four models we have simulated, as illustrated by the subtitles. 



the scalar field Sip in the fundamental observer's frame 
(which is independent of the particle's peculiar velocity 
v) and the part due to frame transformation which is 
proportional to v. In a homogeneous universe, there is 
no spatial gradient and all matter particles are comoving 
with the fundamental observers, so v = 0, and therefore 
the channel II effects vanish. In an inhomogeneous uni- 
verse both terms are nonzero in general, leading to a net 
force which strengthens the attraction between particles, 
and enhances their gravitational clustering. 

The effects arising through channel I appear both 
in homogeneous (via the Friedmann and Raychaudhuri 
equations) and inhomogeneous (via the Poisson equa- 
tion) universes. They arise mainly via a renormalization 
of the contributions to the total density from the cou- 
pled matter species and by a new contribution from the 
scalar field itself. The exact effects depend on the spe- 



cific forms of the scalar potential and coupling function. 
Here, we list our main findings for an inverse power-law 
potential Eq. (13), with an exponential coupling Eq. (14), 
and the model parameters given in Sect. IIIB. 

For the background cosmology, the contribution to the 
Friedmann equation from dark matter is renormalized by 
C(tp) < 1, and C(ip) decreases as the coupling constant 
| -y | increases. As a result, a stronger coupling (larger I7I) 
produces lower cosmic expansion rate during the matter- 
dominated era (Fig. I) 5 , which will in turn helps enhance 
the clustering of matter and promotes structure forma- 



Note in Fig. 1 (lower left panel) the expansion rate of the coupled 
scalar field model finally catches up with that of ACDM at a = 1. 
This is because at late times the dark energy (scalar field) density 
is higher than that in ACDM. 
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FIG. 14: The internal density profiles for dark matter (solid curves) and baryons (dashed curves) separately. We show this for 
the four models we have simulated, as indicated by the subtitles. In each panel the upper two curves are for halo I and the 
lower two curves for halo II. 



tion. 

Such an enhancement is just what we have observed 
from our linear perturbation analysis (Figs. 2 and 3): on 
very large scales, which are beyond the range of the fifth 
force effects, the matter power spectrum increases with 
|7|. On small scales, the fifth force helps to enhance the 
matter power, making it increase further than predicted 
in ACDM (comparing the left to the right panels of Fig. 3, 
where we have separated baryons, which have no scalar 
coupling, from dark matter, which has). This enhance- 
ment starts at a rather early, and so the scalar field docs 
leave imprints on the matter power spectrum at z ~ 50, 
which means that the initial condition for iV-body simu- 
lations also needs to modified. For example, in the model 
with |"y| = 0.15, we have found that : 

(i) The density perturbation at z ~ 49 is about 
10% larger than the corresponding result in ACDM for 



baryons (due to slower cosmic expansion) and about 12% 
larger for dark matter (due to the slower expansion and 
the fifth force effects); 

(ii) The average velocity at the same time is about 
8% larger than the corresponding result in ACDM for 
baryons and about 10% larger for dark matter. We intro- 
duce a quick method to take these changes into account 
when generating initial conditions for iV-body simula- 
tions. 

One of the most important results from our iV-body 
simulations is the magnitude of fifth force. Fig. 6 shows 
that everywhere and at any time (i.e., in both high and 
low density regions) the fifth force (or more precisely, 
the velocity- independent part of it) is proportional to 
gravity in magnitude (with a coefficient 2^ 2 when only 
considering dark matter) . This is the approximation used 
in many previous simplified simulations, and our results 
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confirm numerically that this approximation works fairly 
well, at least for those cases where the scalar potential 
is unimportant, such as ours. Wc emphasize, however, 
that the velocity-dependent part of the fifth force should 
not be dropped (as in some previous simulations) unless 
there is good reason to do so (see [19]). 

For the nonlinear matter power spectrum (Fig. 9), we 
find that the TV-body simulations show agreement with 
linear perturbation analysis on large scales. On interme- 
diate scales, our simulations predict significant enhance- 
ment in the matter power over that predicted both by 
linear perturbation analysis and by the ACDM paradigm. 
This enhancement is strongest at early epochs but grad- 
ually weakens at late times. One possible reason for this 
is that the scalar coupling, C((p), decreases with time, 
and reduces the contribution of dark matter to the grav- 
itational potential, so weakening further clustering, but 
this still needs to be investigated in more detail. 

The bias between dark matter and baryons increases 
with time and with coupling strength \^f\. This is be- 
cause larger I7I implies stronger fifth forces and therefore 
stronger total forces act on dark matter particles com- 
pared to baryons. As time passes, the bias has more time 
to develop and so increases as well (Fig. 10). Using the 
same reasoning, the fifth force increases the attraction of 
dark matter particles and so heavier halos are expected 
to form during the structure formation (Fig. 11). In or- 
der to identify gravitationally bound and virialized halos 
in our models we must also take into account this scalar 
coupling (Sect. IV D). 

The scalar field coupling also has impacts on the in- 
ternal density profiles of the halos (Figs. 12 and 13). We 
average ten of the heaviest halos and find that 

(i) When compared with ACDM, the overdensities in 
the inner regions can be lower while those in the outer 
regions is higher, showing the failure of the halos to retain 
particles in the inner region, either because the particles 
move too fast or because the attractive potential at the 
centre is too weak; 

(ii) The suppression of inner density compared with 
ACDM weakens, rather than strengthens, as I7I increases. 
Indeed, for I7I = 0.2 the density is higher than ACDM 
result throughout the halos. More detailed studies are 
needed to clarify the leading effect responsible for this 
observed pattern. 

(iii) There is a bias between density profiles for baryons 
and dark matter (Fig. 14): the dark matter density is 
always higher because it experiences the fifth force which 
boosts its clustering. 

In summary, in this paper we have been given a com- 
prehensive description of the methodology and imple- 
mentation of general A-body simulations for coupled 
scalar field models. Some important issues in these sim- 
ulations have been addressed here for the first time: the 
consistent solution of the scalar field and fifth force, the 
fifth force effects on generating initial conditions, and the 
effects of scalar field on identifying virialized halos. Al- 
though the situation is complex, wc have identified inter- 



esting new features in these models. We hope these de- 
velopments will lead to more detailed study of nonlinear 
structure formation in this class of models and facilitate 
their subsequent confrontation with observational data. 



Acknowledgments 

The work described in this paper has been performed 
on COSMOS, the UK National Cosmology Supercomputer. 
The coding work at early stage was done on the SARA Su- 
percomputer in the Netherlands, under the HPC-EUR0PA2 
project, with the support of the European Commu- 
nity Research Infrastructure Action under the FP8 Pro- 
gramme. The linear perturbation calculations in this 
work are performed using CAMB [48]. We thank Luca 
Amendola, Marco Baldi, Kazuya Koyama, Andrea Mac- 
cio, Gong-Bo Zhao and Hongsheng Zhao for useful con- 
versations relevant to this work, and computing staff of 
DAMTP for technical help in the usage of COSMOS. B. Li 
is supported by the Research Fellowship at Queens' Col- 
lege, Cambridge, and the STFC. 



Appendix A: Discretized Equations for the TV-body 
Simulations 



In the MLAPM code for Poisson equation Eq. (35) is (and 
in our modified code the scalar field equation of motion 
Eq. (3G) will also be) solved on discretized grid points, 
so wc must develop the discrete versions of Eqs. (33 - 36) 
to be implemented in the code. First, we write down the 
full formalism of the relevant equations. 

Introducing the variable u (cf. Sect. II C), the Poisson 
equation becomes 



V 2 $ r 



— T^CDM i Pc.CDMexp 



+ -f2 B (Pc.B -1) ~ 



7 f Vk<p + 



B 2 Bl 

—'' 

acr 



3Ac 



/-- , B2H S 



3Aa ; 

[y/K0\ 



olV^V 



a, (Al) 



where A is defined in Eq. (20) and is a constant of 0(1) 
(actually A ~ ^de where de means dark energy). Its 
value and that of the quantity \fn,(p are determined solely 
by the background evolution. The computation of the 
background quantities is discussed in Sect. HID. 
The scalar field equation of motion now becomes 



V u = 37fi C DMPc,CDM exp 
3aAa 3 



7 {^J~KLp + 
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3aAa 3 



l + Q 



37^ C DMe 






[Vk<P. 



l+a ' 



(A2) 
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So, in terms of the new variable u, the set of equations 
used in the iV-body code is Eqs. (33, 34) plus Eqs. (Al, 
A2). These equations were used in the code. Among 
them, Eqs. (Al, 34) will use the value of u, while Eq. (A2) 
solves for u. In order that these equations can be inte- 
grated into MLAPM, we need to discretize Eq. (A 2) for the 
application of Newton-Gauss-Seidcl relaxation method. 
This involves writing down a discrete version of this equa- 
tion on a uniform grid with grid spacing h. Suppose we 
want to achieve second-order precision, as is in the de- 
fault Poisson solver of MLAPM, then V 2 m in one dimension 
can be written as 



where a subscript j means that the quantity is evaluated 
on the j'-th point. The generalization to three dimensions 
is straightforward. 

The discrete version of Eq. (A2) is 



L h {ui, jjk ) = 0, 
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in which 
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r 



Then, the Newton-Gauss-Seidel iteration says that we The old solution will be replaced by the new solution to 

can obtain a new (and often more accurate) solution of ^*i,j,fc once the new solution is ready, using the red-black 

u, u°^., using our knowledge about the old (and less Gauss-Seidel sweeping scheme. Note that 
accurate) solution u° ld t via 
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In principle, if we start from a high redshift, then the 
initial guess of u»,j,jfe for the relaxation could be so chosen 
that the initial value of ip in all space is equal to the 
background value (p, because at this time we expect this 
to be approximately true any way. At subsequent time- 
steps we could use the solution for Uj,j,fc from the previous 
time-step as our initial guess. If the timestep is small 
enough then we would expect u to change only slightly 
between consecutive timesteps so that such a guess will 
be good enough for the iterations to converge quickly. 

In practice, however, due to specific features and the 



algorithm of the MLAPM code [39], the above procedure 
may be slightly different in details. 



Appendix B: The Dark Matter Lagrangian 



For simplicity let us write the CDM Lagrangian as 

£cdm = F 1 ^ 2 = V f 9abi a i b (ie., neglecting the mass 
term and ^-function term in Eq. (2)) in which x a = 
dx a /ds and s is some general parameterization of a parti- 
cle's geodesic (we only use this definition in this appendix 
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and note that in other places in the paper x has different 
meaning) . The Euler-Lagrange equation is 



(Bl) 



With these definitions it is straightforward to show that 
Eq. (40) can be expressed as 
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where Ho is the current value of "H, and the coefficients 
of the derivative terms are given by Eqs. (4f , 42) as 
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so that the Lagrange equation becomes 
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For timelike geodesies we can choose s — r (r being the 
proper time) which means that F = 1 so that Eq. (B2) 
becomes the usual geodesic equation 



T bc x x 



= 0. 



(B3) 



If F ^ 1, then we need to retain Eq. (B2) as the geodesic 
equation, which is unfamiliar. We note that F = 1 (and 
£cdm = 1 above) corresponds to the result in our Eq. (4). 
As it stands, Eq. (4) is exact. Since we use proper time r 
for s, the geodesic equation is also written in terms of t, 
which means that when we switch from r to the physical 
time t in Eq. (2,3) there will be some approximation, but 
Eq. (23) is approximate in any case. 



Appendix C: An Algorithm to Solve for the 
Background Evolution 

Here, we aim to give the formulae and an algorithm 
for the background held equations which can also be ap- 
plied to linear Boltzmann codes such as CAMB. Through- 
out this Appendix we use the conformal time r instead 
of the physical time t, and ' = d/dr, % = a' /a. All quan- 
tities appearing here are background ones unless stated 
otherwise. 

For convenience, we will work with dimensionless quan- 
tities and define (p = \fRip and N = In a so that 



$ 



r 



= H 



U 



dip 
dN' 
d 2 ip , dip 



dN 2 



K 



dN' 



(CI) 
(C2) 



In all of the above equations we have used the facts that 

Pcdm oc oT 3 , 

Pb oc oT 3 , 

Prad oc oT 4 . 

We stress again that although coupled to the scalar field, 
the background dark-matter energy density follows the 
same conservation law as in ACDM. The effect of the 
coupling is reflected in the fact that there is a coefficient 
e 7V in front of f^cDM whenever the latter appears in the 
gravitational equations or the scalar field evolution equa- 
tion [17]. 

When solving for p (or <p), we just use Eq. (C3) aided 
by Eqs. (C4, C5). It may appear then that, given any 
initial values for ip- ln i and (dp>/dN) ini , the evolution of tp 
is obtainable. However, Eq. (C4) is not necessarily sat- 
isfied for ip evolved in such way. Instead, it constrains 
the initial condition ip must start with, and the way it 
must subsequently evolve. This in turn is determined by 
the parameters A, a, 7; since a, 7 specify a model and are 
fixed once the model is chosen, the only concern is A. 

For the initial conditions ipi n { and (d<p/dN) ini , we have 
found that the subsequent evolution of (p is rather insen- 
sitive to them. Thus, we choose tp- ln i = (d<p/dN) ini = 
at some very early time (say iV; n i corresponds to a; n ; = 
g-Nini _ io~ 6 ) in all the models. Such a choice is clearly 
not only practical but also reasonable, given the fact that 
we expect that the scalar field starts high up the potential 
and rolls down subsequently. 

As for A, we use a trial-and-error method to find its 
value which ensures that (here a subscript denotes the 
present-day value) 



^rad + ^B 



'CDM 



A 1 f dip 

0% " 1_ 6 [dN 
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which comes from setting a = 1 in Eq. (C4). 

We determine the correct value of A for any given a, 7 
in this way using MAPLE, and then compute the values 
of (p and dip/dN for predefined values of N stored in an 
array. Their values at any time are then obtained us- 
ing interpolation, and with these it is straightforward to 
compute other relevant quantities, such as 'H.,'H! , and if, 
which are used in the Boltzmann and iV-body codes. 



Appendix D: The Zeldovich Approximation 

The initial conditions for iV-body simulations are con- 
ventionally generated using the Zeldovich approximation 
[46, 47], which in its original form works only for non- 
coupled dark matter. When there is a non-minimal cou- 
pling between dark matter particles and a scalar field, it 
must be generalized, and we discuss this here. 

Consider a particle whose actual position r is given as 



a(t)q + b(t)d 



(Dl) 



where q is the Lagrangian coordinate (i.e. initial comov- 
ing coordinate with displacement), d is the displacement 
from q and b(t) governs how the displacement increases 
in time (or, equivalently, how the density perturbation 
grows). From this, we have a deformation tensor 



V: 



dr % 



a(t)5ij + b{t) 



ddj 
d qj ' 



(D2) 



which is just the Jacobian of the coordinate transfor- 
mation, providing information about the change of the 
size of a given volume element (centred on the particle). 
When the deformation is small, as in the case of the early 
times when we set up the initial condition, we have 



det£> 



a 3 (t) 



b 

1 + -V C 

a 



(D3) 



where we have neglected higher-order terms in d, and V q 
is the spatial derivative with respect to the Lagrangian 
coordinate q. Note that the term in brackets is the frac- 
tional change of the volume element due to the collapse 
of matter. 

In a given volume element, because we have two mat- 
ter species now and their density contrasts grow at dif- 
ferent rates because of the different scalar-field coupling, 
it makes sense to have two different 'b(t)' functions: b(t) 
for baryons and b(t) for dark matter. Thus, we end up 
with 



r B = aq + 6d, 
r D = aq + bd, 



(D4) 
(D5) 



in which r B , i"rj denote the actual positions of a baryonic 
particle and a dark-matter particle residing in the volume 
element. We now want to solve for b and b. 



As we will work with constant dark-matter mass, the 
mass conservation is just as simple as in ACDM, and, as- 
suming no particles escape or enter our volume element, 
we have 



p B (t)a 3 d 3 q_ = p B (t,r)d 3 r, 
p D (i)a 3 d 3 q = p D (t,r)d 3 r, 



(D6) 
(D7) 



where p is the average density. Now, from Eq. (D3), it 
follows directly that 
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and so the density contrasts can be expressed in terms of 
b and 6 as 



a 
b. 



-V q d 



-V q d 



■-D+V q -d, 



■-D+Vq-d, 



(D10) 
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where D + , D + will be identified as the linear growth fac- 
tors for baryons and dark matter below. 

Now consider the force laws obeyed by baryons and 
dark matter. For baryons this is very simple. In the New- 
tonian limit: 



r B = -V r 



(D12) 



where V r is the spatial derivative with respect to r, and 
we have 



in the limit of small d. Here, 
tential given by 



aV r (D13) 

4> is the gravitational po- 



: 4ttG [p B + e w po - 2V(f)] , (DI4) 

where we have neglected the radiation, which is negligible 
at later times, and the kinetic energy of the scalar field 
if, which is always small. Meanwhile, the force law for 
dark-matter particles has a contribution from the scalar 
field coupling: 



td = -V r 



7VrV 



(D15) 



where in case that the scalar field potential V(if) is flat 
(as in our case) and the perturbation to the scalar field 
if is small (for redshift z > 50), the scalar field equation 
of motion is approximately 



V^ = 8 7 7rGe w (p D - Pd) 



(D16) 



Writing V = V r from here on, and neglecting the per- 
turbation in if and V(ip), from Eqs. (D10, D14) we have 
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where we have used the Raychaudhuri equation. This can 
be integrated once to obtain 



V4> = r - 4nG (p B D+ + e /v 'p D D+) ad.(D17) 

Eqs. (D4, D12, D17) then give, after some algebra, 



D A 



2-D, 



AttG 



p B D+ + e^poD^ 



0. (D18) 



This is the equation for b. 

The idea of deriving the equation for b is quite similar, 
but now we need to take into account the fifth force. Here, 
Eq. (DIG) can be re-expressed, using Eq. (Dll), as 

VV = SirGa"/e' fV p D D + V ■ d 

which can be integrated once to obtain 

Vip = 87rGa7e w / o D .D+d. (D19) 

Then, Eqs. (D5, D15, D17, D19) lead, again after some 
algebra, to 



Da 



2 a -b + 
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p B D+ + WpbD^ 
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where we have defined j3 = 1 + 2j 2 to encode the effects 
of the fifth force. 

Eqs. (D18, D20) indicate that D +1 D + are the linear 
growth factors of baryons and dark matter. Indeed, if the 
coupling constant 7 = 0, then we find that D + = D+ and 



D, 



2 a -D, 
a 



4nGp m D^ 



0, 



(D21) 



where p m = p B + Pv is the total matter density. This 
is the familiar equation of linear growth equation in this 
approximation. 

In this discussion we have made several approxima- 
tions. For example, the kinetic energy of scalar field, 
which is always sub-dominant, is neglected; the perturba- 
tion of the scalar field potential energy is also neglected 
since the scalar field perturbation is small, especially at 
the time we set up the initial condition; the spatial depen- 
dence of D + is neglected implicitly, because we only con- 
sider small scales where the fifth force simply rescales the 
gravitational constant which governs structure growth. 

Eqs. (D10, Dll) are the starting point of the numer- 
ical codes which generate initial conditions for N-body 
simulations, such as GRAFIC2. Given the matter power 
spectrum at the initial time £,, the code produces density 
fluctuation field <5 as a Gaussian random field and works 
out the displacement field d, or actually -D+d, because 



aq 



fed 



a(q- 



D,d) 



(D22) 



The initial peculiar velocity of the particle is then 



v = ax = aD + d 
where r is the conformal time and 



f ^ D4 



dr 



f = 



dlnD A 

din a 



(D23) 



(D24) 



In our model the initial displacements and velocities of 
particles should be generated separately for the two dif- 
ferent matter species, using their respective matter power 
spectrum and linear growth factor D+ (or D + ). 
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